//////////////////////////////////////////////////////////////////
// (c) Copyright 1998-2002 by Jeongnim Kim
//////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-
#ifndef OHMMS_CP2IS_TEMPLATE_H
#define OHMMS_CP2IS_TEMPLATE_H

/*!\class CP_2I_S 
 * \brief Car-Parrinello dynamics for the electron states
 * Using Tr((2I-S)H) functional for orthonormal basis functions
 */

template<class EigenStates_t, class DensityMatrix_t>
struct CP_2I_S {

  typedef double value_type;

  int cgstep;
  int cgstep_min;
  int cgstep_max;
  int debug_level;
  int cgstep_acc;
  value_type c0, c1, c2;

  EigenStates_t psi0, temp;
  EigenStates_t gr;        //gradient, conjugate gradient

  CP_2I_S() { 
    cgstep = -1; 
    set();
  }

  void 
  set(double em = 134.171, double dt = 1.0, double friction=0.5) {

    // usual verlet if friction == 0, c0 = 2, c1 = -1, c2 = dt*dt/m/2
    double x = dt*friction/2.0;
    c0 = 2.0/(1+x);
    c1 = 1-c0;
    c2 = dt*dt/em/27.2116/(1+x);
  }

  //clean up the vectors
  ~CP_2I_S() { }

  virtual void init(const EigenStates_t& psi) = 0;

  /*!\fn void densitymatrix
   *\brief perform A^t*Q A to caluclate a density matrix
   */
  EigenStates_t&
  Conj(const EigenStates_t& psi) {
    xAB_plus_yC(1.0, Q, psi, 0.0, temp);
    return temp;
  }

  value_type 
  minimize(EigenStates_t& psi, value_type) {

    init(psi); // check the states

    xAtA_plus_yC(-1.0,psi,0.0,Q);
    addDiagonal(Q,2.0);

    apply(psi,temp);

    xABt_plus_yC(1.0, psi, temp, 0.0, H); // H = <Psi|(h|Psi>) = <Psi|Psi_1>

    // G = <Psi|h|Psi> |Psi> - (2I-<Psi|Psi>) h |Psi>
    xAB_plus_yC( 1.0, H, psi, 0.0, gr);
    xAB_plus_yC(-1.0, Q, temp, 1.0, gr);
    
    if(cgstep == 0) {
       psi0 = psi;
       axpy(c2,psi,gr);
    } else {
      copy(psi0,temp);
      copy(psi,psi0);
      scaled(psi,c0);
      axpy(c1,temp,psi);
      axpy(c2,gr,psi);
      //psi = c0*psi0 +c1*temp + c2*gr;
    }
    cgstep++;
    return Trace(H,Q);
  }
};

#endif
/***************************************************************************
 * $RCSfile: CP_2I_S_Solver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: CP_2I_S_Solver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/
